New combined PIC-MCC approach for fast simulation of a radio frequency discharge 

at low gas pressure. 
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A new combined PIC-MCC approach is developed for accurate and fast simulation of a radio 
frequency discharge at low gas pressure and high density of plasma. Test calculations of transition 
between different modes of electron heating in a ccrf discharge in helium and argon show a good 
agreement with experimental data. We demonstrate high efficiency of the combined PIC-MCC 
algorithm, especially for the collisionless regime of electron heating. 

PACS numbers: 52.27. Aj; 52.65.Ww; 52.80.Pi 



I. INTRODUCTION 



The modern trends of plasma technologies are directed 
to a reduction of the gas pressure and an increase of 
plasma density. The further development of efficient 
methods is required for simulation of collisionless regimes 
in a capacitevely coupled and especially in an inductively 
coupled discharges as the collisionless heating of elec- 
trons plays a key role in dynamics of thin skin layers. 
The Particle-in-Cell Monte Carlo Collisions (PIC-MCC) 
method 1] has become a standard simulation technique 
for a gas discharge in plasma reactors of etching or de- 
position. Unlike the fluid approach, the PIC-MCC al- 
gorithm requires larger computer resources, but it pro- 
vides a detail kinetic picture of processes in a gas dis- 
charge. However, a problem of statistical fluctuations 
of an electric field appears at low gas pressures, in par- 
ticular for gases with a deep Ramsauer minimum in the 
elastic scattering cross section. For the periodic electrical 
field E — Eq sin(wt), where lo is the discharge frequency, 
the rate of electron heating in the bulk is proportional to 
vE 2 /{lo 2 + v 2 ), where v is electron collision frequency. At 
high gas pressure in the collisional regime, when v > u>, 
the electric field in the quasineutral plasma is sufficiently 
large and the effect the artificial electron heating is less 
dangerous. At the low gas pressure in the collisionless 
regime (at v < to), the electrons gain the energy in the 
electrode sheaths. In the quasineutral plasma the elec- 
tric field is small and the electrons scattering on the field 
fluctuations essentially distorts the results. Although the 
numerical smoothing of the charge density Q helps to di- 
minish the statistical noise, it is necessary to develop a 
more radical way for reduction of the influence of statis- 
tical fluctuations. 

An interesting idea was suggested in Ref.yj. As the 
discharge simulation lasts more than one thousand of dis- 
charge periods, the averaging of the charge density over 
several periods reduces the statistical noise. But the di- 
rect averaging can lead to the development of the numer- 



ical instability. To eliminate this problem, the electric 
field was calculated in Ref.|3| from the current continu- 
ity equation. However, this approach requires an explicit 
distinction of electrode sheaths, that is difficult for re- 
alization in the two-dimensional case. Besides, it does 
not take into account inertia of electrons, which is very 
important at the low gas pressure. Below we present 
another way of the noise reduction in a new approach 
developed by Vitaly Schweigert. 



II. COMBINED PIC-MCC APPROACH 

In the combined PIC-MCC approach we find the elec- 
tric field distribution from the auxiliary equations which 
are derived from the kinetic equations. The integration 
of the electron and ion kinetic equations over the ve- 
locity gives us the continuity equations for electron and 
ion densities. The integration of the kinetic equations 
multiplied by the velocity gives the continuity equations 
for electron and ion fluxes. The kinetic coefficients are 
calculated with using the electron and ion distribution 
functions, which are found from the electron and ion ki- 
netic equations. To avoid the kinetic coefficients fluctu- 
ations we average them over many periods. Thus, in our 
model the kinetic equations, the auxiliary equations and 
the Poisson equation are solved self-consistently. The ki- 
netic approach allows us to find the kinetic coefficients 
and the electric field distribution is found from the aux- 
iliary equation. The equation system includes the Boltz- 
mann kinetic equations for velocity distribution functions 
of electrons f e (t, x, v) and ions fi(t, x, v), which are three 
dimensional over the velocity and one dimensional in the 
space 
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where v e , Vi, n e , rii, m, M are the electron and ion veloci- 
ties, densities and masses, respectively, E is the electrical 
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field, J e , Ji are the collisional integrals for electrons and 
ions, the transport equations for the density and the flux 
of electrons and ions based on the momentum of the ki- 
netic equations CEJ,© 
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where 
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Q = N g V(Tif e dv e 
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is the ionization rate, a% is the ionization cross sections, 
N g is the gas density, 
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are the effective electron and ion temperatures, respec- 
tively, 

Qe = N g I V ex \v e \<J t f e dVe ~ Ve \ V ex f e dv e , 



Qi = N g J v ix \vi\a r fidvi -i>i J Vi x fidVi 

describe the friction for electrons and ions, the efficient 
frequencies 

_ N g J \v e \a t f e dv e _ N g J \vi\<T r fidVi . . 
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where at is the electron transport cross section, ay is 
the ion resonance charge exchange cross section. Notice 
that in the usual fluid approach the terms Q e , Qi are 
supposed to be zero, which is correct only for the constant 
scattering frequencies. The boundary conditions for the 
auxiliary equations includes the secondary emission as in 
Ref. 4]. It can be easily seen, that the equations -© 
are direct consequences of the kinetic equations GJ,©. 
As we calculate the kinetic coefficients Q, T' e , T/, Q e , 
Qi with solving the kinetic equations JIJ,© with the 
Monte-Carlo method, the obtained densities n' e , n\ have 
to coincide with a good accuracy with values from the 
kinetic equations CQ,@- After calculating the auxiliary 



values of electron n' e and ion densities, we calculate 
the electric field from the Poisson equation 



Ac 
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(10) 



The reduction of statistical noise in our approach is 
reached with averaging the kinetic coefficients Q, T' e , TL 
Q e , Qi over many periods and with smoothing over the 
spatial coordinate. For averaging a function F(x) over 
preceding periods we use the following algorithm 



F(xY = aF{x)' i + (1 - a)F(xY 



(11) 



where F(x) n is the value on the i-period and a = 0.01 
0.1. The spatial smoothing is chosen as in Ref.[2| 



F(x k ) 



F(x k+1 )+2F(x k ) + F(xk- 1 ) 



(12) 



where Xk is the node of the simulation grid. The spa- 
tial smoothing is very important for resolving the space 
charge in the quasineutral part of a discharge, where the 
charge is a small difference of two large and almost equal 
values (ion and electron densities). The computer re- 
sources for solving the transport equations @-@ are 
much smaller than for the kinetic equations, therefore 
the auxiliary equations are solved for each period, and 
the kinetic coefficients are calculated after several peri- 
ods from the kinetic equations Q,J2J- Then the elec- 
tron and ion weights are fitted with the densities n' e , 
n-. We use 5000 simulation particles for each charged 
species, the Cloud-in-Cell charge assignment scheme, the 
null-collisions technique to find the time of electron and 
ion free flight, and the energy conserving scheme with a 
second order of accuracy to solve the equations of mo- 
tion 0,0. The equations ©-© are solved with an 
implicit finite-difference method with using the Schar- 
fetter and Gummel scheme 0. For small grid spacing 
T e » e\cf)k+i — 4>k\ this finite-difference scheme has a sec- 
ond order accuracy in Ax and gives a correct result on 
rough grids for the Boltzmann electron distribution. Like 
for the explicit PIC-MCC method, there exists a restric- 
tion on time step io v A t < 1, where lu p is the plasma 
frequency, for solution of the equations 121-® with the 
Poisson equation |(T0*)| . 

Since the cross section of the electron Coulomb scatter- 
ing is proportional to the electron density and inversely 
proportional to the square of the electron energy, a cor- 
rect discharge simulation of some regimes requires ac- 
counting for electron Coulomb collisions. For description 
of these collisions we apply the method 0, where the 
Langevin force and friction of electrons are introduced 
and defined from their distribution function. Note also, 
that the Coulomb collisions do not change the total elec- 
tron momentum of motion. 

The test calculations show that this algorithm is nu- 
merically stable and allows one to reach a significant 
acceleration of the PIC-MCC method due to two fac- 
tors. At first, the time step for solving the kinetic equa- 
tions with implicit scheme does not depend on 
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the plasma frequency 0, II] ■ At second, averaging over 
many periods allows one to reduce greatly (in 5-20 times) 
the total number of simulation particles in the PIC-MCC 
method without an increase of the statistical noise. 



III. HOW MANY SIMULATION PARTICLES 
WE NEED? 

We study a capacitively coupled radio frequency 
(ccrf) discharge in argon and helium with the combined 
PIC-MCC approach for the experimental conditions of 
Godyak et al\$. We consider an one-dimensional sym- 
metrical ccrf discharge at room temperature for the fre- 
quency v — 13.56 MHz and with the sinusoidal shape 
of the discharge current j. One electrode is grounded 
and the voltage on another electrode is calculated self- 
consistently to provide the desired amplitude of the dis- 
charge current. The spatial grid has typically 81 nodes, 
condensing in electrode sheaths. The minimal grid spac- 
ing is decreased with the gas pressure rise, thus the 
sheath contains approximately the constant number of 
nodes. The cross sections of electron scattering in he- 
lium are taken from ^1 j an d for argon from 0, E3 ■ The 
ion-electron emission on electrodes is taken into account 
with coefficient 0.2 in helium and 0.1 in argon. 

It is known that the statistical error of Monte-Carlo 
methods decreases as 1/N 2 . The statistical noise leads 
to the systematical error in the electron cooling or heat- 
ing. Therefore, first we have studied the impact of the 
number of simulation particles on an accuracy of results 
in three different methods: in the standard PIC-MCC 
0, in the PIC-MCC with the spatial smoothing (PIC- 
MCC SS) @ and in our combined PIC-MCC. The sim- 
ulations are performed for two values of argon pressures 
P = 0.1, 0.3 Torr, the inter-electrode distance d = 2 cm 
and the discharge current j = 2.65 mA/cm 2 . The mean 
electron energy in the discharge center calculated with 
three methods and measured in Ref . |9( is shown in Fig. ^ 
The calculations with the standard PIC-MCC method 
with different N show a significant role of electric field 
fluctuations under the lower gas pressure (squares in 
Fig.^a)). It is seen that the standard PIC-MCC consid- 
erably overestimates the value of e for N = 4000-^256000. 
The second method (PIC-MCC SS) gives much better 
results (circles in Fig. The spatial smoothing indeed 
decreases the statistical noise, but feasibility of this tech- 
nique is restricted, since it distorts the space charge in 
the electrode sheath. At gas pressure P = 0.3 Torr (see, 
Fig. ^b)) when the electron energy e increases with N, 
the PIC-MCC SS method shows the reasonable accuracy 
(within 10%) with small number of simulation particles 
N = 10000. But at the lower pressure P = 0.1 Torr, 
the PIC-MCC SS is not able to provide convergency in 
the electron energy even with N = 256000. It is obvious 
that at low gas pressure in order to obtain the reason- 
able solution with the standard PIC-MCC methods we 
need so an enormous number of the simulation particles 
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FIG. 1: Mean electron energy in the discharge center as 
a function of the total number of simulation particles for 
P = 0.1 Torr (a) and P = 0.3 Torr (b) calculated with the 
standard PIC-MCC method (squares), with the PIC-MCC SS 
method (circles) with spatial smoothing of the space charge 
and electrical field distributions and with our new combined 
algorithm (triangles). 'Cross' is calculation from Ref.Q with 
N = 32000. d = 2 cm, j = 2.65 mA/cm 2 . 
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FIG. 2: Electron density (a) and mean electron energy (b) 
in the discharge center (x = 1 cm) in argon computed (cir- 
cles) and measured in [jj (triangles) for d = 2 cm, j = 
2.65 mA/cm 2 and N = 5000. 
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FIG. 3: Spatial distribution of averaged over period electron 
density (a), potential of electric field (b), mean electron en- 
ergy (c), and electron heating rate (d) in helium for two gas 
pressures P = 0.03 (dashed lines) and 0.3 Torr (solid lines), 
d = 6.7 cm, j = l mA/cm 2 and TV" = 5000. 



FIG. 4: Spatial distribution of averaged over period electron 
density (a), potential of electric field (b), mean electron en- 
ergy (c), and electron heating rate (d) in argon for two gas 
pressures P = 0.03 (dashed lines) and 0.3 Torr (solid lines), 
d = 6.7 cm, j = l mA/cm 2 and N = 5000. 



that these methods are not more applicable. As seen in 
Fig- Hour combined PIC-MCC method gives the electron 
energy which is very close to the experimental one (see, 
Fig. [5J already with the small number of simulation par- 
ticles and the results only weakly depend on N. The 
electron density and the mean electron energy from the 
experiment 0] and from the combined PIC MCC simu- 
lations with N = 5000 are shown in Fig. [21 as a function 
of gas pressure. The calculated dependence of the e from 
P demonstrates the transition between different modes 
of the electron heating found in and well agrees with 
the experimental data. 

IV. VALIDITY OF THE COMBINED PIC-MCC 

APPROACH. SIMULATION RESULTS OF A 
CCRF DISCHARGE IN HELIUM AND ARGON. 

Depending on gas pressure, there are two different 
regimes of electron heating (collision and collisionless) 
in rf discharges which are well studied experimentally 
and numerically (see, for example [lL l9l fill IT3 . fl3. Il4| ^ . 



The collision electron heating takes place due to elastic 
scattering of the electrons on the atoms, when the di- 
rected velocity transfers into the thermal one. At high 
gas pressures the collisonal (or ohmic) heating controls 
the electron energy in the quasineutral part of the dis- 
charge. At the low gas pressure the electrons are heated 
due to interaction with moving sheaths boundaries and 
the ohmic heating in bulk is very small. For these two 
regimes the spatial distributions of the electron density, 
the electrical potential, the mean electron energy and the 
electron heating rate Wh = —eE J v ex f e dv e are shown in 
Figs. in helium and argon, respectively. 

The results are obtained for two different gas pressures 
P = 0.03 Torr and P = 0.3 Torr, for d — 6.7 cm and 
j = 1 mA/cm 2 . As expected, in helium the mean elec- 
tron energy increases with pressure lowering in order to 
compensate an increase of particle losses at the electrodes 
and in argon we observed the opposite behavior. The rea- 
sons of reduction of the electron e nerg y under pressure 
lowering in argon are discussed in [la |. where a drop of 
the electron temperature up to the gas temperature is 
predicted in the absence of Coulomb electron collisions. 
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FIG. 5: Electron energy probability function in helium in the 
discharge center (x = 3.35 cm) for different gas pressures, 
d = 6.7 cm, j = 1 mA/cm 2 and N = 5000. 



(b) 




e (eV) 

FIG. 6: Electron energy probability function in argon in the 
discharge center (x = 3.35 cm) for different gas pressures, 
d = 6.7 cm, j = 1 mA/cm 2 and N = 5000. 

Note, that in helium larger heating rate (in the center 
of discharge) refers to lower e. This non-local effect can 
not be predicted within the fluid or the diffusion-drift 
approaches. The electron energy probability functions 
(EEPF) are shown in Figs. 15161 for helium and argon. 
The data presented in these figures averaged over the dis- 
charge period. As in the experiment we also found in 
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FIG. 7: Electron-neutral elastic cross sections in argon (solid 
line) and in helium (dotted line) as functions of the electron 
energy. 



argon that the EEPF changes from a Druyvesteyn shape 
to a bi-Maxwellian one with decreasing the gas pressure 
(see, Fig. IHJ. At the low gas pressure the electrons are 
separated into two groups. The cold electrons are not 
able to reach the sheath boundary and their ohmic heat- 
ing is very weak due to Ramsauer minimum in the elas- 
tic cross section (see, Fig. 0. The fast electrons heated 
in the sheaths maintain the discharge operation and pro- 
vide the gas ionization. Fig.[5]presents the computed and 
measured j9( electron temperature (T e = 2U e /3) in the 
discharge center (x = 3.35 cm). The decrease of the gas 
pressure is accompanied with a drop of e. A comparison 
with experimental data shows a good agreement (within 
20 -j- 30%) within a pressure range P — 0.03 -j- 0.3 Torr 
for helium and for argon. 

The calculation gives higher energy at the larger gas 
pressure P = 1 Torr. The difference between computed 
and measured data at higher gas pressure is likely due 
to the contribution of metastable states in the ionization 
kinetics, especially in helium (see, for example 01 )• ^ n 
the model of electron-neutral collisions in our simulations 
we do not take into account the multi-step ionization. At 
low gas pressures we have better agreement because the 
metastable atoms are deactivated on electrodes and the 
influence of multi-step ionization reduces. The study of 
ionization kinetics in noble gases is out of the scope of 
this article. Note that in our earlier study of the ccrf 
discharge in helium we have considered the metastable 
atoms and obtained a good agreement with experimental 
data for high gas pressures. 

In conclusion we have presented the combined PIC- 
MCC approach for fast simulation of the rf discharge 
over a wide range of gas pressures and current densities. 
The validity of the new approach is justified by compar- 
ison with the experiment data. The advantage of our 
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approach is the considerable decrease of the number of 
simulation particles N. We are able to reach a speed- 
up factor of ten for the collision regime and even more 
for the collisionless regime compared with the standard 
PIC-MCC calculations. 



FIG. 8: Effective electron temperature (T e = 2U e /3) in the 
discharge center (x — 3.35 cm) in helium (a) and argon (b). 
Computed T e (circles) and measured T e 9] (triangles) for d = 
6.7 cm, j = 1 A/cm 2 and N = 5000. 
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